# CLEAN WORKSPACE AND LOAD PACKAGES --------------------------------------------

rm(list = ls())
library(datasim)
library(tidyverse)

# SIMULATE MULTIVARIATE SPATIAL DATA -------------------------------------------

# set.seed(4)
Corr <- matrix(c(1, -0.3, 0, -0.3, 1, 0.3, 0, 0.3, 1), nrow = 3)
sigmas <- rep(0.4^0.5, 3)
D <- diag(sigmas)
Cov <- D %*% Corr %*% D

# beta <- c(-0.5, 0, 0.5)
beta <- c(0, 0, 0)
variance <- 0.6 * matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1), nrow = 3)
cor.model <- "exp_cor"
cor.params <- list(list(phi = 0.05), list(phi = 0.05), list(phi = 0.1))

f <- list(
  mean ~ mfe(x1, beta = get("beta")) +
    mre(factor(id), sigma = get("Cov")) +
    mgp(list(s1), variance = get("variance"), cor.model = get("cor.model"),
        cor.params = get("cor.params"), range = 2),
  sd ~ I(0)
  )

n <- 300
m <- 3
(data_geo <- sim_model(formula = f, n = n, responses = m))
## # A tibble: 900 x 9
##       id     x1    s1 mre.factor.mean mgp.list.mean   mean    sd response
##    <int>  <dbl> <dbl>           <dbl>         <dbl>  <dbl> <dbl>    <dbl>
##  1     1  0.791 0.733          0.772        -0.311   0.461     0    0.461
##  2     2  0.659 1.53          -0.374        -0.432  -0.806     0   -0.806
##  3     3 -0.514 1.16          -0.697        -1.06   -1.75      0   -1.75 
##  4     4 -0.826 1.34           0.285        -0.677  -0.392     0   -0.392
##  5     5 -1.05  1.01           0.522         0.135   0.657     0    0.657
##  6     6 -0.636 0.689          0.653        -0.267   0.386     0    0.386
##  7     7 -0.283 1.73          -0.0330       -0.806  -0.839     0   -0.839
##  8     8 -0.636 0.239          0.124        -0.363  -0.239     0   -0.239
##  9     9  0.141 0.964          0.113        -0.331  -0.218     0   -0.218
## 10    10  0.682 0.670         -0.432         0.0897 -0.342     0   -0.342
## # ... with 890 more rows, and 1 more variable: response_label <int>
# knitr::kable(head(data_model, 10))

# VISUALIZE MULTIVARIATE SPATIAL DATA ------------------------------------------

ggplot(data_geo, aes(x1, response)) +
  geom_smooth(aes(col = factor(response_label))) +
  geom_point(aes(col = factor(response_label)))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

ggplot(data_geo, aes(s1, mgp.list.mean)) +
  geom_line(aes(col = factor(response_label)))

vg <- data_geo %>%
  mutate(s2 = s1) %>%
  group_by(response_label) %>%
  nest() %>%
  mutate(variog = purrr::map(data, ~ gstat::variogram(mgp.list.mean ~ 1, ~ s1 + s2, . ,
                                                      cutoff = 3, width = 0.005))) %>%
  dplyr::select(-data) %>%
  unnest()

ggplot(vg, aes(dist, gamma)) +
  geom_point(aes(size = np)) +
  geom_smooth() +
  expand_limits(y = 0, x = 0) +
  scale_x_continuous(limits = c(0, 3)) +
  scale_y_continuous(limits = c(0, 1.3)) +
  stat_function(fun = function(x) 0.6 * (1-exp(- x/0.05)), col = 2, size = 1) +
  facet_wrap(~response_label, ncol = 1)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 85 rows containing non-finite values (stat_smooth).
## Warning: Removed 85 rows containing missing values (geom_point).

vg <- data_geo %>%
  mutate(s2 = s1) %>%
  group_by(response_label) %>%
  nest() %>%
  mutate(variog = purrr::map(data, ~ gstat::variogram(response ~ 1, ~ s1 + s2, . ,
                                                      cutoff = 3, width = 0.005))) %>%
  dplyr::select(-data) %>%
  unnest()

ggplot(vg, aes(dist, gamma)) +
  geom_point(aes(size = np)) +
  geom_smooth() +
  expand_limits(y = 0, x = 0) +
  scale_x_continuous(limits = c(0, 3)) +
  scale_y_continuous(limits = c(0, 2)) +
  stat_function(fun = function(x) 0.4 + 0.6 * (1-exp(- x/0.05)), col = 2, size = 1) +
  stat_function(fun = function(x) 0.4 + 0.6 * (1-exp(- x/0.1)), col = 3, size = 1) +
  facet_wrap(~response_label, ncol = 1)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 48 rows containing non-finite values (stat_smooth).
## Warning: Removed 48 rows containing missing values (geom_point).

# vg <- data_geo %>%
#   mutate(s2 = s1) %>%
#   dplyr::filter(response_label == 3) %>%
#   gstat::variogram(mgp.list.mean ~ 1, ~ s1 + s2, . , cutoff = 3, width = 0.005)
# ggplot(vg, aes(dist, gamma)) +
#   geom_point(aes(size = np)) +
#   geom_smooth() +
#   expand_limits(y = 0, x = 0) +
#   scale_x_continuous(limits = c(0, 3)) +
#   stat_function(fun = function(x) 0.6 * (1-exp(- x/0.05)), col = 2, size = 2)


data_geo %>%
  dplyr::select(id, mre.factor.mean, response_label) %>%
  spread(response_label, mre.factor.mean) %>%
  dplyr::select(-id) %>%
  GGally::ggpairs(aes(fill = "any"))

data_geo_wide <- data_geo %>%
  dplyr::rename(ability = response, id_person = id) %>%
  gather(var, value, mre.factor.mean:ability) %>%
  mutate(var = paste0(var, response_label)) %>%
  select(-response_label) %>%
  spread(var, value)

# SIMULATE ITEM FACTOR DATA ----------------------------------------------------

q <- 10
init_data <- purrr::map(1:q, ~ data_geo_wide) %>%
  purrr::reduce(rbind)

# n <- 300
difficulty <- matrix((1:q - 5)/10 * 2, nrow = 1)
discrimination1 <- seq(0.4, 1.5, length.out = q)
discrimination2 <- runif(q, 0, 2)
discrimination3 <- runif(q, 0, 2)
discrimination1[1] <- 1
discrimination1[c(3, 5, 8)] <- 0
discrimination2[1:2] <- c(0, 1)
discrimination2[c(4, 5, 10)] <- 0
# discrimination3[1:3] <- c(0, 0, 1)
# discrimination1 <- discrimination1 * 0.3
# discrimination2 <- discrimination2 * 0.3
cbind(discrimination1, discrimination2, discrimination3)
##       discrimination1 discrimination2 discrimination3
##  [1,]       1.0000000        0.000000      0.40052452
##  [2,]       0.5222222        1.000000      1.10412733
##  [3,]       0.0000000        1.098587      1.32626150
##  [4,]       0.7666667        0.000000      1.57900160
##  [5,]       0.0000000        0.000000      1.08222874
##  [6,]       1.0111111        1.783167      1.32947141
##  [7,]       1.1333333        1.505002      1.74563860
##  [8,]       0.0000000        0.865549      0.05820813
##  [9,]       1.3777778        0.300353      0.52441984
## [10,]       1.5000000        0.000000      0.70927425
f <- list(
  prob ~ mfa(ones, beta = get("difficulty")) +
    mfe(ability1, beta = get("discrimination1")) +
    mfe(ability2, beta = get("discrimination2")),
  # + mfe(ability3, beta = get("discrimination3")),
  size ~ I(1)
  )

data_long <- sim_model(formula = f,
                        link_inv = list(pnorm, identity),
                        generator = rbinom,
                        responses = q,
                        n = n,
                        init_data = init_data
                        )

data_long <- dplyr::rename(data_long, subject = id,
                           item = response_label, y = response)

# VISUALIZE ITEM FACTOR DATA ---------------------------------------------------

explor <- data_long %>%
  group_by(subject) %>%
  summarize(endorse = mean(y),
            ability1 = unique(ability1),
            ability2 = unique(ability2),
            # ability3 = unique(ability3),
            x1 = unique(x1))
ggplot(explor, aes(ability1, endorse)) + geom_point(alpha = 0.5)

ggplot(explor, aes(ability2, endorse)) + geom_point(alpha = 0.5)

# ggplot(explor, aes(ability3, endorse)) + geom_point(alpha = 0.5)
# ggplot(explor, aes(x1, endorse)) + geom_point(alpha = 0.5)

# PREPARE DATA -----------------------------------------------------------------

response <- data_long$y
coordinates <- dplyr::select(data_geo_wide, s1)
dist <- as.matrix(dist(coordinates))
# dist <- as.matrix(dist(dplyr::select(data_geo_wide, s1)[order(data_geo_wide$s1),]))
# dist <- dist[order(data_geo_wide$s1),]
n
## [1] 300
q
## [1] 10
m <- 2
iter <- 5 * 10 ^ 4
thin <- 5
# iter <- 5 * 10 ^ 3
cor.params <- c(0.04, 0.04)
sig.params <- c(0.6 ^ 0.5, 0.6 ^ 0.5)
fix.sigma <- 0.4^0.5
# sigma_prop <- matrix(c(0.138, -0.023, -0.023, 0.1), 2) * 2.38 ^ 2 / 2
sigma_prop <- 0.001 * diag(5)
disc_mat <- cbind(discrimination1, discrimination2)
L_a <- lower.tri(disc_mat, diag = TRUE) * 1
L_a[c(3,5,8), 1] <- 0
L_a[c(4,5,10), 2] <- 0
T_gp <- diag(m)
# diag(T_gp) <- 0
# T_gp[2,2] <- 0

# RUN --------------------------------------------------------------------------

Rcpp::sourceCpp("../src/mirt-gibss-sp.cpp")
source("../R/ggplot-mcmc.R")
Rcpp::sourceCpp("../src/ifa-driver.cpp")
source("../R/spmirt.R")

# # set.seed(5)
# system.time(
#   samples <- ifa_gibbs_sp(response, dist, n, q, m, cor.params, sig.params,
#                           Corr[1:2, 1:2], fix.sigma, sigma_prop, L_a, T_gp, 0.234,
#                           iter)
# )

iter <- 5000
thin <- 5
system.time(
  samples <- spmirt(
    response = response,  predictors = NULL, coordinates = coordinates,
    standardize = TRUE,
    nobs = n, nitems = q, nfactors = 2, niter = iter, thin = thin,
    constrains = list(A = L_a, W = T_gp, V_sd = sigmas[1:2]/2),
    adaptive = list(Sigma = NULL, Sigma_R = NULL, Sigma_gp_sd = NULL,
                    Sigma_gp_phi = NULL, scale = 1, C = 0.7, alpha = 0.8,
                    accep_prob = 0.234),
    sigmas_gp_opt = list(initial = 0.6, prior_mean = 0.6, prior_sd = 0.4),
    phi_gp_opt = list(initial = 0.05, prior_mean = 0.05, prior_sd = 0.4))
  )
## Standardixing
##    user  system elapsed 
## 437.705 305.496 205.306
iter = iter / thin
thin2 <- 1

attr(samples, "model_info")[-c(1, 2, 3)]
## $nobs
## [1] 300
## 
## $nitems
## [1] 10
## 
## $nfactors
## [1] 2
## 
## $ngp
## [1] 2
## 
## $niter
## [1] 5000
## 
## $thin
## [1] 5
## 
## $standardize
## [1] TRUE
## 
## $constrain_L
##       [,1] [,2]
##  [1,]    1    0
##  [2,]    1    1
##  [3,]    0    1
##  [4,]    1    0
##  [5,]    0    0
##  [6,]    1    1
##  [7,]    1    1
##  [8,]    0    1
##  [9,]    1    1
## [10,]    1    0
## 
## $constrain_T
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
## 
## $constrain_V_sd
## [1] 0.3162278 0.3162278
## 
## $adap_Sigma
##       [,1]  [,2]  [,3]  [,4]  [,5]
## [1,] 0.001 0.000 0.000 0.000 0.000
## [2,] 0.000 0.001 0.000 0.000 0.000
## [3,] 0.000 0.000 0.001 0.000 0.000
## [4,] 0.000 0.000 0.000 0.001 0.000
## [5,] 0.000 0.000 0.000 0.000 0.001
## 
## $adap_scale
## [1] 1
## 
## $adap_C
## [1] 0.7
## 
## $adap_alpha
## [1] 0.8
## 
## $adap_accep_prob
## [1] 0.234
## 
## $c_initial
##  [1]  0.6984590 -0.8918420 -0.3100677  0.3747685  0.5509894 -0.1521639
##  [7] -0.2630168  0.1391695 -0.4260893 -0.4334037
## 
## $c_prior_mean
##  [1] 0 0 0 0 0 0 0 0 0 0
## 
## $c_prior_sd
##  [1] 1 1 1 1 1 1 1 1 1 1
## 
## $A_initial
##       [,1] [,2]
##  [1,]    1    0
##  [2,]    0    1
##  [3,]    0    0
##  [4,]    0    0
##  [5,]    0    0
##  [6,]    0    0
##  [7,]    0    0
##  [8,]    0    0
##  [9,]    0    0
## [10,]    0    0
## 
## $A_prior_mean
##       [,1] [,2]
##  [1,]    1    0
##  [2,]    0    1
##  [3,]    0    0
##  [4,]    0    0
##  [5,]    0    0
##  [6,]    0    0
##  [7,]    0    0
##  [8,]    0    0
##  [9,]    0    0
## [10,]    0    0
## 
## $A_prior_sd
##       [,1] [,2]
##  [1,] 0.45 1.00
##  [2,] 1.00 0.45
##  [3,] 1.00 1.00
##  [4,] 1.00 1.00
##  [5,] 1.00 1.00
##  [6,] 1.00 1.00
##  [7,] 1.00 1.00
##  [8,] 1.00 1.00
##  [9,] 1.00 1.00
## [10,] 1.00 1.00
## 
## $R_initial
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
## 
## $R_prior_eta
## [1] 1.5
## 
## $B_initial
##      [,1] [,2]
## [1,]   NA   NA
## 
## $B_prior_mean
##      [,1] [,2]
## [1,]   NA   NA
## 
## $B_prior_sd
##      [,1] [,2]
## [1,]   NA   NA
## 
## $sigmas_gp_initial
## [1] 0.6 0.6
## 
## $sigmas_gp_mean
## [1] 0.6 0.6
## 
## $sigmas_gp_sd
## [1] 0.4 0.4
## 
## $phi_gp_initial
## [1] 0.05 0.05
## 
## $phi_gp_mean
## [1] 0.05 0.05
## 
## $phi_gp_sd
## [1] 0.4 0.4
## 
## $model_type
## [1] "spifa"
samples_tib <- as_tibble(samples, iter/2)
#summary(samples_tib)
samples_long <- gather(samples_tib)

as_tibble.spmirt.list(samples, 0, thin2, "c") %>%
  gg_trace(alpha = 0.6)

as_tibble.spmirt.list(samples, 0, thin2, "a") %>%
  gg_trace(alpha = 0.6)

as_tibble.spmirt.list(samples, iter/2, thin2, "a") %>%
  gg_density(alpha = 0.5, ridges = TRUE, aes(fill = Parameters), scale = 4)
## Picking joint bandwidth of 0.122

as_tibble.spmirt.list(samples, iter/2, thin2, "theta") %>%
  dplyr::select(1:100) %>%
  gg_density(alpha = 0.5, ridges = TRUE, aes(fill = Parameters), scale = 4)
## Picking joint bandwidth of 0.13

as_tibble.spmirt.list(samples, 0, thin2, "theta") %>%
  select(1:10) %>%
  gg_trace(alpha = 0.6)

as_tibble.spmirt.list(samples, 0, thin2, "corr") %>%
  gg_trace(alpha = 0.6)

as_tibble.spmirt.list(samples, 0, thin2, "mgp_sd") %>%
  gg_trace(alpha = 0.6)

as_tibble.spmirt.list(samples, iter/2, thin2, "mgp_sd") %>%
  gg_density(alpha = 0.5, ridges = FALSE, aes(fill = Parameters), scale = 4) +
  stat_function(fun = dlnorm, colour = "red",
                args = list(meanlog = log(0.6), sdlog = 0.4))
## Warning: Ignoring unknown parameters: scale

as_tibble.spmirt.list(samples, 0, thin2, "mgp_phi") %>%
  gg_trace(alpha = 0.6)

as_tibble.spmirt.list(samples, iter/2, thin2, "mgp_phi") %>%
  gg_density(alpha = 0.5, ridges = FALSE, aes(fill = Parameters), scale = 4) +
  stat_function(fun = dlnorm, colour = "red",
                args = list(meanlog = log(0.05), sdlog = 0.4))
## Warning: Ignoring unknown parameters: scale

as_tibble.spmirt.list(samples, 0, thin2, "a") %>%
  gg_density2d(`Discrimination 1`, `Discrimination 2`, each = 10,
               keys = c("Item ", "Discrimination "),
               highlight = c(discrimination1, discrimination2))
## Warning: Computation failed in `stat_density2d()`:
## bandwidths must be strictly positive
## Warning: Computation failed in `stat_density2d()`:
## bandwidths must be strictly positive

## Warning: Computation failed in `stat_density2d()`:
## bandwidths must be strictly positive

## Warning: Computation failed in `stat_density2d()`:
## bandwidths must be strictly positive

## Warning: Computation failed in `stat_density2d()`:
## bandwidths must be strictly positive

## Warning: Computation failed in `stat_density2d()`:
## bandwidths must be strictly positive

as_tibble.spmirt.list(samples, 0, thin2, "a") %>%
  gg_scatter(`Discrimination 1`, `Discrimination 2`, each = 10,
               keys = c("Item ", "Discrimination "),
               highlight = c(discrimination1, discrimination2))

as_tibble.spmirt.list(samples, iter/ 2, select = "a") %>%
  summary() %>%
  mutate(param = c(discrimination1, discrimination2)) %>%
  gg_errorbarh() +
  geom_point(aes(param, Parameters), col = 3)

as_tibble.spmirt.list(samples, iter/2, select = "c") %>%
  summary() %>%
  mutate(param = as.numeric(difficulty)) %>%
  gg_errorbarh() +
  geom_point(aes(param, Parameters), col = 3)

as_tibble.spmirt.list(samples, iter/2, select = "theta") %>%
  dplyr::select(1:300) %>%
  summary() %>%
  mutate(param = data_geo$response[1:300]) %>%
  gg_errorbarh(sorted = TRUE) +
  geom_point(aes(x = param), col = 3)

as_tibble.spmirt.list(samples, iter/2, select = "theta") %>%
  dplyr::select(301:600) %>%
  summary() %>%
  mutate(param = data_geo$response[301:600]) %>%
  gg_errorbarh(sorted = TRUE) +
  geom_point(aes(x = param), col = 3)

ability1_pred <- as_tibble.spmirt.list(samples, iter/2, select = "theta") %>%
  dplyr::select(1:300) %>%
  summary() %>%
  mutate(param = data_geo$response[1:300],
         s1 = data_geo$s1[1:300],
         s2 = s1,
         estim = `50%`)
ability1_pred %>%
    ggplot(aes(s1, `50%`)) +
    geom_line() +
    geom_line(aes(s1, param, col = "real"))

vg <- gstat::variogram(estim ~ 1, ~ s1 + s2, ability1_pred, cutoff = 3, width = 0.01)
ggplot(vg, aes(dist, gamma)) +
  geom_point(aes(size = np)) +
  geom_smooth() +
  expand_limits(y = 0, x = 0) +
  scale_x_continuous(limits = c(0, 3))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

ability1_pred <- as_tibble.spmirt.list(samples, iter/2, select = "theta") %>%
  dplyr::select(1:300)
ability1_pred <- ability1_pred[nrow(ability1_pred),]
ability1_pred <- ability1_pred %>%
  summary.spmirt() %>%
  mutate(param = data_geo$response[1:300],
         s1 = data_geo$s1[1:300],
         s2 = s1,
         estim = `50%`)
ability1_pred %>%
    ggplot(aes(s1, `50%`)) +
    geom_line() +
    geom_line(aes(s1, param, col = "real"))

vg <- gstat::variogram(estim ~ 1, ~ s1 + s2, ability1_pred, cutoff = 3, width = 0.01)
ggplot(vg, aes(dist, gamma)) +
  geom_point(aes(size = np)) +
  geom_smooth() +
  expand_limits(y = 0, x = 0) +
  scale_x_continuous(limits = c(0, 3))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

ability2_pred <- as_tibble.spmirt.list(samples, iter/2, select = "theta") %>%
  dplyr::select(301:600) %>%
  summary() %>%
  mutate(param = data_geo$response[301:600],
         s1 = data_geo$s1[301:600],
         s2 = s1,
         estim = `50%`)
ability2_pred %>%
  ggplot(aes(s1, `50%`)) +
  geom_line() +
  geom_line(aes(s1, param, col = "real"))

vg <- gstat::variogram(estim ~ 1, ~ s1 + s2, ability2_pred, cutoff = 3, width = 0.005)
ggplot(vg, aes(dist, gamma)) +
  geom_point(aes(size = np)) +
  geom_smooth() +
  expand_limits(y = 0, x = 0) +
  scale_x_continuous(limits = c(0, 2))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 160 rows containing non-finite values (stat_smooth).
## Warning: Removed 160 rows containing missing values (geom_point).

ability2_pred <- as_tibble.spmirt.list(samples, iter/2, select = "theta") %>%
  dplyr::select(301:600)
ability2_pred <- ability2_pred[nrow(ability2_pred),]
ability2_pred <- ability2_pred %>%
  summary.spmirt() %>%
  mutate(param = data_geo$response[301:600],
         s1 = data_geo$s1[301:600],
         s2 = s1,
         estim = `50%`)
ability2_pred %>%
  ggplot(aes(s1, `50%`)) +
  geom_line() +
  geom_line(aes(s1, param, col = "real"))

vg <- gstat::variogram(estim ~ 1, ~ s1 + s2, ability2_pred, cutoff = 3, width = 0.005)
ggplot(vg, aes(dist, gamma)) +
  geom_point(aes(size = np)) +
  geom_smooth() +
  expand_limits(y = 0, x = 0) +
  scale_x_continuous(limits = c(0, 2))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 160 rows containing non-finite values (stat_smooth).

## Warning: Removed 160 rows containing missing values (geom_point).

# # # PREPARE DATA FOR MODELLING ---------------------------------------------------
# #
# # Y <- data_model %>% dplyr::select(id, response, response_label) %>%
# #   spread(response_label, response) %>%
# #   arrange(id) %>%
# #   dplyr::select(-id) %>%
# #   as.matrix()
# #
# # X <- data_model %>% dplyr::select(id, matches("^x[[:digit:]]+$")) %>%
# #   unique() %>%
# #   arrange(id) %>%
# #   dplyr::select(-id) %>%
# #   as.matrix()
# #
# # Beta <- matrix(beta, nrow = 1)
# # Sigma_proposal <- diag(1, 3)
# #
# # # RUN MODEL --------------------------------------------------------------------
# #
# # getwd()
# # Rcpp::sourceCpp("../src/multi-lm.cpp")
# # source("../R/ggplot-mcmc.R")
# #
# # iter <- 10^6
# # system.time(
# #   samples <- multi_lm(Y, X, iter, 0.01 * Sigma_proposal, 0.001 * Sigma_proposal)
# # )
# # samples %>% map(~ tail(.))
# #
# # # apply(samples$beta, 2, mean)
# # # cor(samples$beta)
# #
# # # Visualize traces
# # as_tibble(samples, 0, 100, select = "beta") %>%
# #   gg_trace(wrap = TRUE, alpha = 0.6)
# #
# # as_tibble(samples, 0, 100, select = "beta") %>% gg_trace(alpha = 0.6)
# # as_tibble(samples, 0, 100, select = "corr_chol") %>% gg_trace(alpha = 0.6)
# # as_tibble(samples, 0, 100, select = "corr") %>% gg_trace(alpha = 0.6)
# # as_tibble(samples, 0, 100, select = "sigmas") %>% gg_trace(alpha = 0.6)
# #
# # bla <- as_tibble(samples, iter/2, select = "sigmas")
# # cov(log(bla))
# # nrow(unique(bla)) / nrow(bla)
# #
# # bla <- as_tibble(samples, iter/2, select = "corr_chol")
# # cov(bla)
# # nrow(unique(bla)) / nrow(bla)
# #
# # # Visualize densities
# #
# # as_tibble(samples, iter / 2, select = "corr_chol") %>%
# #   gg_density(aes(fill = Parameters), scale = 2, alpha = 0.5, ridges = TRUE)
# #
# # as_tibble(samples, iter / 2, select = "corr") %>%
# #   gg_density(aes(fill = Parameters), scale = 1, alpha = 0.5, ridges = TRUE)
# #
# # # Visualize credible intervals
# # as_tibble(samples, iter / 2, select = "beta") %>%
# #   summary() %>%
# #   mutate(param = beta) %>%
# #   gg_errorbarh() +
# #   geom_point(aes(param, Parameters), col = 3)
# #
# # Corr_chol <- t(chol(Corr))
# # corr_chol <- Corr_chol[lower.tri(Corr_chol, diag = TRUE)]
# # corr <- Corr[lower.tri(Corr)]
# #
# # as_tibble(samples, iter / 2, select = "corr_chol") %>%
# #   summary() %>%
# #   mutate(param = corr_chol) %>%
# #   gg_errorbarh() +
# #   geom_point(aes(param, Parameters), col = 3)
# #
# # as_tibble(samples, iter / 2, select = "corr") %>%
# #   summary() %>%
# #   mutate(param = corr) %>%
# #   gg_errorbarh() +
# #   geom_point(aes(param, Parameters), col = 3)
# #
# #
# # as_tibble(samples, iter / 2 ,select = "sigmas") %>%
# #   summary() %>%
# #   mutate(param = sigmas) %>%
# #   gg_errorbarh() +
# #   geom_point(aes(param, Parameters), col = 3)
# #
# #
# # # Visualize credible intervals for all Parameters
# # as_tibble(samples, iter / 2) %>%
# #   summary() %>%
# #   mutate(param = c(beta, corr_chol, corr, sigmas)) %>%
# #   gg_errorbar() +
# #   geom_point(aes(Parameters, param), col = 3)
# #